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1 Introduction 

This chapter will very briefly introduce and review some computational 
experiments in using trainable gene regulation network models to simulate and 
understand selected episodes in the development of the fruit fly. Drosophila 
melanogaster. For details the reader is referred to the papers introduced below. It 
will then introduce a new gene regulation network model which can describe 
promoter-level substructure in gene regulation. 

As described in chapter 2, gene regulation may be thought of as a combination of 
cis-acting regulation by the extended promoter of a gene (including all 



regulatory sequences) by way of the transcription complex, and of trans-acting 
regulation by the transcription factor products of other genes. If we simplify the 
cis-action by using a phenomenological model which can be tuned to data, such 
as a unit or other small portion of an artificial neural network, then the full trans- 
acting interaction between multiple genes during development can be modelled 
as a larger network which can again be tuned or trained to data. The larger 
network will in general need to have recurrent (feedback) connections since at 
least some real gene regulation networks do. This is the basic modeling 
approach taken in (Mjolsness et al. 1991), which describes how a set of recurrent 
neural networks can be used as a modeling language for multiple developmental 
processes including gene regulation within a single cell, cell-cell communication, 
and cell division. Such network models have been called "gene circuits", "gene 
regulation networks", or "genetic regulatory networks", sometimes without 
distinguishing the models from the actual modeled systems. 

In (Mjolsness et al. 1991) a number of choices were made in formulating the 
trainable gene regulation network models, which affect the spatial and temporal 
scales at which the models are likely to be useful. The dynamics was chosen to 
operate deterministically and continuously in time, on continuous-valued 
concentration-like variables, so that the dynamical equations for the network are 
coupled systems of ordinary differential equations (ODE's). One such form was 



dv f } 

in which v, is the continuous-valued state variable for gene product i, T tj is the 

matrix of positive, zero, or negative connections by which one transcription 
factor can enhance or repress another, and g() is a nonlinear monotonic sigmoidal 
activation function. When a particular matrix entry T i} is nonzero, there is a 

regulatory "connection" from gene product j to gene i . The regulation is 
enhancing if T is positive and repressing if it is negative. If 7 is zero there is no 

connection. Figure 1 sketches the model, drawing a few representative nonzero 
connections as arrows between "genes" represented by open circles. The entire 
network is localized to a cell but communicates with other such networks in 
nearby cells. 

[Figure 1 goes here.] 

Such equations are often stiff due to the nonlinear transfer function g(u). 
Optimizing the unknown parameters T, and _ has so far proven to be 
computationally difficult: special versions of simulated annealing optimization 
(Lam and Delosme 1988a, 1988b) have been required for good results, e.g. to start 
from expression patterns derived from a known model and recover its 
parameters reliably (Reinitz and Sharp 1995). As discussed in Chapter 2, this 
kind of training is quite different and much slower than the usual 



"backpropation of error" training used with feed-forward (nonrecurrent) 
artificial neural networks. Informatics work on improving this situation could be 
important. 

In addition to the analog circuit model, the framework of (Mjolsness et al. 1991) 
also proposes a dynamic grammar by which multiple biological mechanisms can 
be modeled by networks and then combined into a consistent overall dynamical 
model. The grammar/ circuit combination has some similarities to hybrid 
models incorporating Discrete Event Systems and ODE's. In this way one can for 
example combine intracellular and intercellular regulation network submodels. 
The grammar is also suitable for implementation in object-oriented computer 
simulation programs. 


2. Three Case Studies 

In this section, some of the literature on trainable gene circuit models which have 
been fit to Drosophila gene expression patterns is reviewed. Three applications to 
pattern formation are shown to demonstrate the generality of the methods. First, 
a model of gap gene expression patterns along the anterior-posterior axis will be 
described. Second, the extension of this model to incorporate the important pair- 
rule gene eve and a number of other improvements will be introduced. Finally, a 


gene circuit model of neurogenesis incorporating nonlinear signaling between 
nearby cells through the Notch receptor and the Delta ligand will be briefly 
described. 



2.1 Gap Gene Expression 


Such gene regulation network models can be tuned or "trained" with real gene 
expression data, and then used to make robust and at least qualitatively correct 
experimental predictions, as was shown in (Reinitz et al. 1992). In that study the 
goal was to understand the network of gap genes expressed in bands (domains) 
along the anterior-posterior (A-P) axis of the very early embryo (the syncytial 
blastoderm) of Drosophila. This experimental system has the advantage that there 
are no cell membranes between adjacent cell nuclei, so elaborate cell-cell 
signalling mechanisms do not need to be modeled. Also Drosophila is an easy 
species to manipulate genetically, as for example "saturation mutagenesis" - 
finding all the genes affecting a particular process - is possible. 

Positional information along the A-P axis of the syncytial blastoderm is encoded 
in a succession of different ways during development. At first the main 
encoding is a roughly exponential gradient of bicoid (bed) protein imposed by the 
mother fly, along with maternal hunchback (hb) expression. These provide gene 
regulation network inputs to the gap genes: Kruppel (Kr), knirps (kni), giant (gt), 
tailless (til), and hunchback (hb) again. These each establish one or two broad 
domains of expression along the A-P axis. The gap genes then serve as network 
inputs to the pair-rule genes including even-skipped (eve) and fushi tarazu (ftz ), 
which establish narrow, precise stripes of expression and precise positional 



coding. These in turn provide input to segment-polarity genes such as engrailed 
and wingless which are the first to retain their expression pattern into adulthood. 
For example, engrailed is expressed in bands just one cell wide which define the 
anterior borders of the parasegments. Introductions to the relevant Drosophila 
developmental biology may be found in (Lawrence 1992) and (Alberts et al. 
1994). 

An example of a spatial gene expression pattern along the A-P axis of a triple- 
stained embryo is shown in Figure 2. Here, fluorescently labelled antibodies 
simultaneously label those nuclei in the syncytial blastoderm expressing Kruppel, 
giant, and even-skipped. 


[Figure 2 goes here.] 

The first computer experiments with fitting such analog gene regulation nets to 
real expression data concerned the establishment of the broad gap gene domains 
(excluding the extreme ends of the A-P axis) from maternally supplied initial 
conditions, by a gene regulation network in which all gap genes interact with all 
others and bed provides input to, but does not receive any input from, the gap 


genes. 


Figure 3 shows the experimentally observed and model-fitted curves for gap 
gene expression. They are in qualitative agreement, which is the most that can 
be expected from the expression data that was available at the time. The extra 
dip in gt expression could not be predicted by the model, which can be 
interpreted as an indication of the role of circuit components not included in the 
model. 


[Figure 3 goes here.] 

The most important predictions of the model concerned the anomalous dose- 
response observed by (Driever and Nusslein-Volhard 1988). Figure 4 shows the 
prediction in detail; it may be summarized by saying that positional information 
for the gap gene system is specified cooperatively by maternal bed and hb. This 
qualitative behavior was observed to be robust over many runs of the simulated 
annealing parameter-fitting procedure, and therefore taken to be a prediction of 
the model. Essential features of the cooperative control of positional information 
by maternal bed and hb were verified experimentally in (Simpson-Brose et al. 
'1994). The gap gene model prediction and the experiment ocurred 
independently of one another. 


[Figure 4 goes here.] 


2.2 Eve Stripe Expression 


Following the gap gene computer experiments, (Reinitz and Sharp 1995) went on 
to perform a detailed study of the gap gene circuit as extended to include the 
first of the pair-rule genes, eve. The further observations which could be 
included in this model allowed an important milestone to be reached: not only 
qualitative behaviors, but also the circuit parameter signs and rough magnitudes 
became reproducible from one optimization run to another, and some 
parameters such as connections to eve were still more reproducible. Hence, far 
more could be predicted. For example the diffusion constant for eve was much 
lower than for other transcription factors in successful runs. This has an 
experimental interpretation: eve mRNA is expressed in the outer part of each 
future cell just as the cell membranes are invaginating into the blastoderm 
embryo, providing an apical obstruction to diffusion. 

More importantly, each of the eight boundaries of the four central stripes of eve 
expression could be assigned a particular gap gene as the essential controller of 
that boundary. This picture is in agreement with experimental results with the 
possible exception of the posterior border of eve stripe 3, the interpretation of 
which is an interesting point of disagreement (Small et al. 1996, Reinitz and 



Sharp 1995, Frasch and Levine 1987) and a possible focal point for further 
laboratory and/ or computer experiments. 

Further experimental understanding of the gap genes' influence on eve 
expression is obtained in (Reinitz et al. 1998), where it is shown that the fact that 
eve is unregulated by other pair-rule genes can be understood by the phase of its 
periodic spatial pattern: no other phase offset pattern of pair-rule expression (e.g. 
the phase-shifted patterns of hairy or fushi-tarazu) can be produced from gap gene 
input alone. 

[Figure 5 goes here.] 

Related work on modeling the gap gene and eve system of A-P axis positional 
information in Drosophila includes (Hamahashi and Kitano 1998). 


2.3 Neurogenesis and Cell-Cell Signaling 

The syncytial blastoderm is very favorable, but also very unusual, as 
morphogenetic systems go because there is no cell membrane interposed 
between nearby cell nuclei and therefore the elaborate mechanisms of cell-cell 
signaling do not come into play. But if we are to model development in its 
generality it is essential to include signaling along with gene regulation 



networks. As a first attempt in this direction, we have modeled the selection of 
particular cells in an epithelial sheet (later in Drosophila development) to become 
neuroblasts. Virtually the same gene network is thought to be involved in the 
selection of particular cells in wing imaginal disks to be sensor organ precursors. 
The essential molecule to add is the Notch receptor, a membrane-bound receptor 
protein responsible for receiving the intercellular signals which mediate this 
selection process. It binds to a ligand molecule ("Delta" for this system) on 
neighboring cells. Recent experiments (Schroeter et al. 1998) indicate that it acts 
on the nucleus (following activation by a ligand on another cell) by having an 
intracellular domain cleaved off and transported there. Variants of the Notch 
receptor occur in many developmental subsystems where a subpopulation of 
cells must be picked out, in Drosophila and homologously across many species. 

In (Marnellos 1997) and (Marnellos and Mjolsness 1998a, 1998b) are reported 
computer experiments incorporating both intracellular and intercellular 
components in a gene regulation network model of neurogenesis. A minimal 
gene circuit model with lateral inhibition (such as depicted in Figure 6) was not 
quite sufficient to produce the observed patterns of selection robustly. 
Incorporating a denser intracellular connection matrix and/or the dynamic 
effects of delamination on the geometry of cell/ cell contact area produced better 
results. However, the "data" to which the fits were made was highly abstracted 
from real gene expression data so it is premature to draw a unique biological 



hypothesis from the model. Figure 7 shows the resulting model behavior in the 
case of dense interconnections. 

[Figure 6 goes here.] 

Related work on Notch-mediated signaling in Drosophila developmental models 
includes the appearance of Notch and Delta in the ommatidia model of 
(Morohashi and Kitano 1998). 

[Figure 7 goes here.] 


3 Extending the Modeling Framework to Include Promoter Substructure 

A very important scientific problem is to understand the influence of promoter 
substructure on eve stripe formation. The eve promoter has many transcription 
factor binding sites, some of which are grouped more or less tightly into 
promoter elements such as the stripe 2 "minimal stripe element" (MSE 2) (Small 
et al. 1992), or a similar less tightly clustered element for stripes 3 and 7 (Small et 
al. 1996). As an example of the scientific problems that are raised, if is hb an 
enhancer for MSE 2 but an inhibitor for MSE 3, what is its net effect on eve and 
can it change sign (Reinitz et al. 1998)? And how are we to understand the action 
of "silencer" elements such as the one apparently responsible for long-range 



repression of zen by dorsal (Gray et al. 1995)? Such questions point to the need 
for at least one additional level of complication in the phenomenological models 
of gene networks whose application is described above, to describe the 
substructure of promoters: binding sites, their interactions, and promoter 
elements. Otherwise the relevant experiments cannot even be described, let 
alone predicted, with network models. 

In (Small et al. 1992) an informal model for activation of MSE 2 is suggested: it is 
activated by bed and hb "in concert", and repressed by gt anteriorly and Kr 
posteriorly. A simple "analog logic" expression for the activation of MSE 2 in 
terms of variables taking values in [0,1] might then be (GRN 1998): 

u mse 2 = { bcd + 7 x hb)( 1 - £t)(l - Kr) 

V MSE2 = S{ U MSEl) 

where _ is a weight on the relative contribution of hb vs bed. A similar simplified 
formula for the model of (Small et al. 1996, figure 8) for MSE 3 could be for 
example: 


u mse 2 = Dstat(\ - hb){ 1 - kni ) 

V MSEJ = S{ U MSE2 ) 



(We omit direct activation of MSE3 by tailless (til) since til represses kni (Pankratz 
et al. 1989) which represses MSE3.) The rate of eve transcription would be 
approximated by a further analog logic formula including a weighted "or" of the 
MSE activations v MSE2 and v MJ£3 . 

The validation or invalidation of such formulae and their interpretation in terms 
of more detailed models will require a quantitative treatement of the relevant 
expression data which is not yet available. It may also lead to fitting the 
parameters in quantitative network models of promoter-level substructure 
within a gene regulation network. 


3.1 An Example: Hierarchical Cooperative Activation 

As an example of such a gene network model incorporating promoter level 
substructure, I introduce here a "Hierarchical Cooperative Activation" (HCA) 
model for the degree of activation of a transcription complex. It at least seems 
more descriptive of known mechanisms than a previous attempt to derive 
phenomenological recurrent neural network equations as an approximation to 
gene regulation dynamics (Mjolsness et al. 1991). An earlier suggestion for 
including promoter-level substructure in gene regulation networks is described 
in (Sharp et al. 1993). The present HCA model is more detailed but has not been 



fit to any experimental data yet and is therefore quite speculative: perhaps a next 
stage of successful modeling will include some of the following ingredients. 

The basic idea of the model is to use an equilibrium statistical mechanics model 
(complete with partition functions valid for dilute solutions (Hill 1985)) of 
"cooperative activation" in activating a protein complex. Such a model can be 
constructed from the following partition function, which is essentially the 
Monod-Wyman-Changeux (MWC) model for a concerted state change among 
subunits (Hill 1995): 


z=*n<i+v /lw >+n<>+v, w ) 

b b 

in which the probability of activation of some complex is determined by relative 
binding constants for each component b of the complex in the active and inactive 
states, but there are no other interactions. As before, v. represents the 

concentration of gene product j of a gene circuit. In this formula, j is a function 
of b so that each binding site is specialized to receive only one particular 
transcription factor. To remove this assumption one could write instead 


b j b i 

where A = 0 or 1 specifies which transcription factors may bind to which sites by 
its sparse nonzero elements. For either expression, K is the relevant binding 



constant for a binding site when the complex is in its "active" state and K is the 
binding constant when the complex is inactive. 

For this partition function, given a global active or inactive state, all binding sites 
are independent of one another. For example the components could be the 
occupants of all the binding sites b within a particular regulatory region of a 
eukaryotic promoter. This conditional independence leads to the products over 
the binding sites in the expression for Z. There are two such products because 
there is one additional bit of global state which can be "active" or "inactive". 

For this model the probability of activation of the complex under consideration 
can be calculated and it is: 


P = 8(u) = 


Ku 

1 + Ku 


»-n 


1 + Kb v m ) 


so (if Kv « 1) 


u 




(Further simplifications result if the binding constants K b specific for a given 
transcription factor j(b) are all roughly equal to a common value K i . The final 
line above suggests a neural-network like approximation for u, although in that 



regime g could be linearized also.) We will use this model as a building block to 
construct a more detailed one. 

Given the MWC-style model of "cooperative activation", we'd like to use it 
hierarchically: to describe the activation of promoter "modules" or "elements" in 
terms of transcription factor concentrations, and then again to describe the 
activation of the whole transcription complex in terms of the "concentrations" of 
active promoter elements, which are proportional to their activities. An 
additional wrinkle is to allow either monomeric or homodimeric transcription 
factor binding. (Heterodimers will be introduced with appropriate notation 
later.) The resulting bare-bones hierarchical model would replace the neural-net 
activation dynamics 


T, — L = [transcribing]; - A ( v, 
dt 

[transcribingl = g(u,) 
u i= y L T IJ v j +h i 


j 


with the two-level model 


T i — = [transcribing]; - A f -v, 
dt 

Ju { 


[transcribing^ = g(u ( ) = y— ^ 


Ju ; 


“,=n 




\ +7 p 

aei V 1 ~ J a I a 



(the product is taken over enhancer elements which regulate genei) and 


p = o (u ) = ^2— 

1 + K a u a 


n 

bea 


f \ + v v" (ft) ^ 

1 + A fc V ;(fr) 

1 4 - Y 1 f n ^ 

1 + K b v i(h) J 


Here n(b)-l for monomers and 2 for homodimers. Note that for this simple feed- 
forward version of the model, the parameters K b and K h are related to 
observables 


f ah = 


V \j n ^ 

i + v;,» 


fah = 


Y \) n ^ 

i+v«> 


where is the probability that site b e is occupied if _ is active, and is the 

probability that site b is occupied if _ is inactive. In principle these quantities 
could be observed by in vivo footprinting. Such observations could be used to 
evaluate the parameters K h to use in the first expression for u a for arbitrary 
inputs Vj. 

If we are modeling a network rather than a single gene, then some of the 
quantities listed above require an additional i index. 



We have the opportunity to include a few more important biological mechanisms 
at this point. One is the possibility that, as in the Endol6 model of (Yuh et al. 98), 
the hierarchy could go much deeper than two levels - especially if transcription 
complex formation is a sequential process. Another significant mechanism is 
competitive binding within a promoter element. This could arise if several 
transcription factors bind to a single site, as we have formulated earlier, or if 
binding at one site eliminates the possibility of binding at a nearby site and vice 
versa. In this case the 4-term product of two 2-term binding-site partition 
functions is replaced with one three-term function by excluding the 
configuration in which both competing sites are occupied: 

z»- = 0 

j k 

4- = (i + Z aA v / + 

i * 

(where again A = 0 or 1 describes which transcription factors bind to which sites 
by its sparse nonzero elements) with corresponding modifications to the update 
equations. Also homodimeric and heterodimeric transcription factor binding are 
easy to accommodate with appropriate concentration products in more general 
one-site and two-site partition functions: 



7 (l) . 

= d 

+ XWy 

1 

+ X\t^ v y v *) 
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= 0 

+ I^Vi 
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Transcription factor trimers and higher order subcomplexes at adjacent binding 
sites could be described by suitable generalizations of these expressions, at the 
cost of introducing more parameters. 


Similarly, constitutive transcription factor binding with activation by 
phosphorylation or dephosphorylation can be described with minor 
modifications of the appropriate one-site or two-site partition functions. For 
example one could use Michaelis-Menton kinetics in steady-state for 
phosphorylation and dephosphorylation, and the one-site dimeric partition 
functions would become 


4" = a + x.» 

* = (1 + ^ A bjk , lm ^bfk V j V k X ! l( X l + 3m)) 
y* 


where x, is proportional to the concentration of a kinase for the bound j/k dimer 
(with proportionality constants depending on the catalytic reaction rates) and y m 



is proportional to a corresponding phosphatase concentration. Also A hJk lm < A hjk , 

so that the extra indices / and m just specify the relevant kinase(s) and 
phosphatase(s) from a kinase network. For example MAP kinase mediated 
signaling could be modeled as activating a gene regulation network by this 
mechanism. 


In this model formulation we have omitted lateral interactions other than 
competitive binding between activation of nearby binding sites. Such 
interactions could be modeled in the manner of an Ising model. For simplicity 
we just use a tree topology of states and partition functions here. 


Given such one-site and two-site partition functions, the overall partition 
function for a promoter element in terms of its binding sites is: 


Z a = K a 


n? 

V*|C=0 
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Here each binding site competes with at most one other one as determined by the 
0/1-valued parameters C b ,C w . 


In this picture, silencers are just particular promoter elements with sufficiently 
strong negative regulation of transcription to veto any other elements. 



The Hierarchical Cooperative Activation (HCA) dynamics then become 


T, — = [ transcribing ], - A ( v, 
dt 

[transcribing], = gU) = Ju ‘ 

1 + JU: 
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with Z's as before: 
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and likewise for inactive-module (hatted) Z's and K's. These partition functions 
encode monomeric, homodimeric and heterodimeric protein-DNA binding using 
the various A parameters. 



The resulting HCA model (Figure 8) can describe promoter elements, silencer 
regions, dimeric and competitive binding, and constitutive transcription factor 
binding, among other mechanisms. The price is that there are considerably more 
unknown parameters in the model than in the previous recurrent neural network 
models - not exponentially many as in the general N-binding site partition 
function, but enough to pose a challenge to model-fitting procedures and data 
sets. 


4 Conclusion 

Gene regulation networks have been applied to model several episodes in the 
development of Drosophila, successfully making contact with experimental 
results. A variety of biological mechanisms including intercellular signaling can 
now be included in such models. We proposed a new version of gene regulation 
network models for use in describing experiments which involve promoter 
substructure, such as transcription factor binding sites or promoter regulatory 


elements. 
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Figure 1. Sketch of recurrent analog neural network model for gene regulation 
networks. A set of analog-valued units v are connected in a continuous-time, 
recurrent circuit by a connection matrix T. Communication with circuits in other 
cells may require additional connections, e.g. as formulated in [Mjolsness et al 
'91]. 

Figure 2. Spatial pattern of gene expression in a Drosophila syncytial blastoderm 
for two gap genes and one pair-rule gene. Immunofluorescent staining of nuclei 
for Kruppel (green), giant (blue), and even-skipped (red). Overlap areas of Kr and 
eve appear yellow, and overlaps of gt and eve appear purple. Image courtesy of 
John Reinitz. 

Figure 3. Data and model for gap gene circuit. Horizontal axes are nuclei along 
lateral midline from anterior to posterior. Vertical axes are relative 
concentrations. Left: data estimated from immunofluorescence images similar to 
Figure 2 for pairs of gap genes. Right: output of a circuit model fit to expression 
data using a nonlinear least squares criterion and simulated annealing 
optimization. From [Reinitz et al. '92], 

Figure 4. Predictions of the model as bicoid dosage is increased: location of 
selected landmarks along A-P (horizontal) axis vs. number of bed copies (vertical 
axis), (a) Displacement of a landmark (anterior margin of the Kr domain) 



expected if it were determined by reading off a fixed concentration value of 
maternal Bicoid protein alone, (b-c) Smaller displacement of the same landmark 
(anterior margin of the Kr domain) predicted by model. (A retrodiction.) (d) 
Observed anomalously small displacement of a related landmark: the first eve 
stripe, not available in the gap gene model but expected to be offset anteriorly 
from the Kr landmark. Note anomalously high slope compared to a, but as in 
b,c. (e) Prediction: return to the behavior of (a) if maternal hunchback is set equal 
to zero. From [Reinitz et al. '92], 

Figure 5. Drosophila eve stripe expression in model (right) and data (left). 

Green: eve expression, red: kni expression. From [Reinitz and Sharp '95], 
Courtesy J. Reinitz and D. H. Sharp. 

Figure 6. A hypothesized minimal gene regulation circuit for lateral inhibition 
mediated by Notch and Delta. Redrawn from [Heitzler et al '96, Figure 6]. Two 
neighboring cells express Notch (N) and Delta (DI) at their surfaces. Notch 
positively regulates transcription of genes of the Enhancer-of-split complex 
E(spl)-C, which negatively regulate transcription of genes of the achaete-scute 
complex (AS-C), which positively regulate transcription of Delta. Curved 
boundaries are the cell membranes between two neighboring cells. Related 
circuit diagrams have been suggested elsewhere e.g. [Lewis '96]. 



Figure 7. Cluster resolution. A circuit "trained" to resolve simple proneural 
cluster configurations into individual neuroblasts (or sensory organ precursor 
cells) is tested on more complex and irregular configurations. In this case each 
cluster was successfully resolved into a single neuroblast, but the large clusters 
resolve more slowly. Times: t=l (top left), t=76 (top right), t=106 (bottom left), 
t=476 (bottom right). Similar to [Mamellos and Mjolsness '98a]; courtesy George 
Mamellos. 

Figure 8. Hierarchical Cooperative Activation (HCA) model for promoter 
substructure within a gene "node" in a gene regulation network. Different layers 
of sub-nodes have different forms of dynamics. This network could be used to 
selectively expand some or all of the nodes in Figure 1, for example just the "eve" 
gene in a network for the gap genes and eve. 
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